rm(list = ls())

setwd("/Users/Joan/Google Drive/Article Twitter Streams and Removal Rates/Data")
load("complete_data205k.RData")
library(psych)
library(plyr)
library(arm)
library(lme4)
library(Hmisc)

m1 <- glmer(removed ~ (1|.id), family=binomial (link="logit"), data=all.tweets)
summary(m1)
m2 <- glmer(removed2 ~ (1|.id), family=binomial (link="logit"), data=all.tweets)
summary(m2)
for(i in 1:length(coef(m2)[[1]][,1])){
  coefs <- coef(m2)[[1]][i,]
  curve(as.double(coefs[1])*x, 
        from=-20, to=20, col=i, add=TRUE)
  
}
## MODEL 1: Plot random intercepts

u0 <- ranef(m1, condVar=T)
u0se <- sqrt(attr(u0[[1]], "postVar")[1,,])
commid <- as.numeric(as.factor(rownames(u0[[1]])))
u0tab <- cbind("commid" = commid, "u0" = u0[[1]], "u0se" = u0se)
colnames(u0tab)[2] <- "u0"
u0tab <- u0tab[order(u0tab$u0),]
u0tab <- cbind(u0tab, c(1:dim(u0tab)[1]))
u0tab <- u0tab[order(u0tab$commid),]
colnames(u0tab)[4] <- "u0rank"
u0tab <- u0tab[order(u0tab$u0rank, decreasing = F),]
names <- rownames(u0tab)

u1 <- ranef(m2, condVar=T)
u1se <- sqrt(attr(u1[[1]], "postVar")[1,,])
commid <- as.numeric(as.factor(rownames(u1[[1]])))
u1tab <- cbind("commid" = commid, "u1" = u1[[1]], "u1se" = u1se)
colnames(u1tab)[2] <- "u1"
u1tab <- u1tab[order(u1tab$u1),]
u1tab <- cbind(u1tab, c(1:dim(u1tab)[1]))
u1tab <- u1tab[order(u1tab$commid),]
colnames(u1tab)[4] <- "u1rank"
u1tab <- u1tab[order(u1tab$u1rank, decreasing = F),]
names2 <- rownames(u1tab)

tiff("Fig5.tif", width=2300, height=2300, pointsize=12, res=300)

par(mfrow=c(1,2), oma = c(0,2,2,0), mar=c(3,5,3,1)) #B, L, U, R
plot(u0tab$u0, u0tab$u0rank, type="n", axes=F, frame.plot=F,
     main="Removed Before 15 min", ylab='', font.main=1, cex.main = 0.8,
     xlab="",xlim=c(-1.8, 1.8))
segments(u0tab$u0 - 1.96*u0tab$u0se, u0tab$u0rank, u0tab$u0 + 
           1.96*u0tab$u0se, u0tab$u0rank)
points(u0tab$u0, u0tab$u0rank, cex=0.3, col="blue")
abline(v=0, col="red")
axis(1, at=c(-1.8,0, 1.8), label=c("-1.8", "0", "1.8"))
axis(2, at=(1:36), labels=c(names) , par(las=1))
text(u0tab$u0, u0tab$u0rank, round(u0tab$u0, digits = 3), offset = 0.2, cex=0.4, pos=3, col="blue")

plot(u1tab$u1, u1tab$u1rank, type="n", axes=F, frame.plot=F, 
     main="Removed Before 7 days", ylab='', font.main=1, cex.main = 0.8,
     xlab="", xlim=c(-2.5, 2.5))
segments(u1tab$u1 - 1.96*u1tab$u1se, u1tab$u1rank, u1tab$u1 + 
           1.96*u1tab$u1se, u1tab$u1rank)
points(u1tab$u1, u1tab$u1rank, cex=0.3, col="blue")
abline(v=0, col="red")
axis(1, at=c(-2.5,0, 2.5), label=c("-2.5", "0", "2.5"))
axis(2, at=(1:36), labels=c(names2) , par(las=1))
text(u1tab$u1, u1tab$u1rank, round(u1tab$u1, digits = 3), offset = 0.2, cex=0.4, pos=3, col="blue")

mtext("Random Intercepts By Term, Null Model",outer = TRUE, font=3, cex = 1.2)

dev.off()

